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ABSTRACT 



Using the results of an extended set of high-resolution non-radiative hydrodynamic sim- 
, ulations of galaxy clusters we obtain simple analytic formulae for the dark matter and hot gas 

J> ■ distribution, in the spherical approximation. Starting from the dark matter phase-space radial 

l/"") | density distribution, we derive fits for the dark matter density, velocity dispersion and veloc- 

. ity anisotropy. We use these models to test the dynamical equilibrium hypothesis through the 

Jeans equation: we find that this is satisfied to good accuracy by our simulated clusters in- 
side their virial radii. This result also show that our fits constitute a self-consistent dynamical 
model for these systems. 

We then extend our analysis to the hot gas component, obtaining analytic fits for the gas 
density, temperature and velocity structure, with no further hypothesis on the gas dynamical 
status or state equation. Gas and dark matter show similar density profiles down to w 0.06R V 
Q-f (with R v the virial radius), while at smaller radii the gas flattens producing a central core. 

q Gas temperatures are almost isothermal out to roughly Q.2R V , then steeply decrease reaching 

;_i ■ at the virial radius a value almost a factor of two lower. We find that the gas is not at rest 

^ \ inside R v : velocity dispersions are increasing functions of the radius, motions are isotropic 

CS • to slightly tangential, and contribute non-negligibly to the total pressure support. We test this 

model using a generalization of the hydrostatic equilibrium equation, where the gas motion 
is properly taken into account. Again we find that the fits provide an accurate description of 
the system: the hot gas is in equilibrium and is a good tracer the overall cluster potential if all 
terms (density, temperature and velocity) are taken into account, while simpler assumptions 
cause systematic mass underestimates. In particular, we find that using the so-called p-model 
underestimates the true cluster mass by up to 50 per cent at large radii. We also find that, if 
gas velocities are neglected, then a simple isothermal model fares better at large radii than a 
non-isothermal one. 

The shape of the gas density profile at small radii is at least partially explained by the 
gas expansion caused by energy transfer from dark matter during the collapse. In fact, when 
gas bulk energy is also considered, gas and dark matter are in energy equipartition in the final 
system at radii r > 0.1R V , while at smaller radii the gas is hotter than the dark matter. This 
energy unbalance is also probably the reason of the further global halo compression compared 
to a pure collisionless collapse, which we point out by comparing the dark matter and total 
density profiles of our hydro-simulated clusters with a set of identical - but pure N-body - ones. 
The compression has the effect of rising the mean concentration by an amount of roughly 10 
per cent. 

Key words: cosmology: theory - galaxies: clusters - hydrodynamics - X-rays: galaxies - 
dark matter - methods: numerical 



1 INTRODUCTION 

Galaxy clusters are the largest virialized systems we observe in the 
universe today. To first order they can be described as dark matter 
(DM) haloes, since dark matter makes up at least 80 per cent of 
* rasia@pd.astro.it their total mass. The second most important component of galaxy 
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clusters is a hot plasma, called the intracluster medium (ICM); it 
has temperature of the order of a few keV, roughly tracing the depth 
of the gravitational potential wells of the systems. 

Any comparison with observations in the optical or via grav- 
itational lensing requires modelling at least the dark matter com- 
ponent of galaxy clusters, while X-ray or millimetric observations 
[e.g. the Sunyaev-Zel'dovich (SZ) effect] require knowledge of the 
internal structure of the ICM. Both aspects can in principle be stud- 
ied by mean of numerical simulations. In fact, in the last few years 
simulations have reached enough resolution to reliably describe 
the structural properties of DM haloes down to the kpc scale. An 
important result along this line is the discovery that relaxed dark 
matter haloes exhibit a density profile which can be accurately 
described by a double power-law, with inner asymptotic slope of 
around —1 and outer asymptotic slope close to —3 (Navarro, Frenk 
& White 1996, 1997, hereafter NFW; Cole & Lacey 1996; Tormen, 
Bouchet & White 1997). Although the exact value of the exponent 
of the inner slope is still somewhat debated (Moore et al. 1998; Jing 
& Suto 2000; Jing 2000; Ghigna et al. 2000; Power et al. 2003), the 
overall model is quite well defined and robust. Fitting formulae for 
the density profile of galaxy clusters are very useful to probe the 
standard paradigm of structure formation and cosmological mod- 
els against observations, like the distribution of galaxies in clusters 
(e.g. the CNOC project: Carlberg et al. 1997), and the gravitational 
lensing properties of galaxy clusters (e.g. Bartelmann et al. 1998; 
Meneghetti, Bartelmann & Moscardini 2003). 

However, density profiles contain information only on the 
structural properties of dark matter haloes, while knowledge of the 
velocity structure is required if one wants to extend the descrip- 
tion to the dynamical properties. A few attempts have been done 
in this direction, both using analytical calculations and numeri- 
cal simulations (e.g. Sheth et al. 2001; Taylor & Navarro 2001; 
Hiotelis 2002a,b). These papers usually made simplifying assump- 
tions on the dynamics (e.g. isotropic velocity dispersion or a given 
functional form for the phase-space density or for the velocity 
anisotropy) and imposed the Jeans equation as a constraint to derive 
the model equations. Still, the difficulty of the problem has made 
it impossible to generate so far a self-consistent analytical model 
of the density and velocity structure of dark matter haloes with no 
assumption on the density or velocity field besides spherical sym- 
metry. An easier - but equally interesting - way is that of taking 
the answer directly from numerical simulations, and use the mean 
radial profiles of a sample of simulated clusters as a first guess for 
such dynamical model, with no assumption on the dynamical status 
of the system. This is the first goal of the present paper. 

The second, and perhaps even more interesting goal, is extend- 
ing this attempt to modelling the hot intracluster gas. Given a dou- 
ble power-law functional form for the dark matter density profile, it 
is possible to calculate the radial distribution of the hot gas by im- 
posing the hydrostatic equilibrium equation and assuming that the 
gas is either isothermal or polytropic (e.g. Suto, Sasaki & Makino 
1998; Komatsu & Seljak 2001; Ascasibar et al. 2003, hereafter 
AYMG03). Very recently, Lee & Suto (2003) presented a physical 
model for the non-spherical distribution of the intracluster gas in 
hydrostatic equilibrium under the gravity of triaxial dark matter ha- 
los (Jing & Suto 2002): their solutions for the gas and temperature 
profiles have been obtained in both isothermal and polytropic equa- 
tions of state. In the present work we instead take the approach of 
deriving all relevant quantities directly from the simulations, with 
no assumptions on the gas state equation or dynamical equilibrium. 
We measure the mean radial gas density and temperature profiles 
of our simulated clusters and model them using analytical fits. As 



the gas in simulated clusters is not completely at rest inside the 
virial radius, we also derive fits for the gas velocity dispersions and 
anisotropy. The fits are then tested against the hydrostatic equilib- 
rium equation - and not imposing it. The resulting model gives us 
an idea of how relaxed are, on average, simulated clusters and how 
well they are described by the equilibrium equations. We believe 
that it is important to have a complete model for ICM+DM because 
observations measure the total mass distribution, and the presence 
of gas, even if not dominant, can well affect the DM component. 

The plan of the paper is as follows. In Section 2 we present the 
hydrodynamical N-body simulations used in this work. In Section 
3 we model the dark matter density and the velocity structure, and 
then use our model to test the dynamical equilibrium of clusters 
through the Jeans equation, cast in the form of a mass estimator. In 
Section 4 we perform a similar analysis for the hot gas component. 
In Section 5 we estimate the effect of the presence of ICM on the 
density profile. Moreover we discuss the dependence of our results 
on the dynamical and environmental status of galaxy clusters and 
we compare our results to previous work. Finally, our conclusions 
are summarized in Section 6. 



2 NUMERICAL SIMULATIONS 

In order to produce our sample of simulated galaxy clusters, we 
used the re-simulation technique, based on a software package de- 
veloped by our group (ZIC: Zoom Initial Conditions). The method 
is presented in Tormen et al. (1997). Here we give only a brief sum- 
mary of the adopted procedure. 

We started from a cosmological N-body simulation (the 'par- 
ent simulation'), with 512 3 particles in a box of side 479/z~' Mpc. 
This has been produced by N. Yoshida for the Virgo Consortium 
(Yoshida, Sheth & Diaferio 2001; see also Jenkins et al. 2001). 
The assumed cosmological model is a flat universe, where the con- 
tributions to the density parameter from dark matter, baryons and 
cosmological constant are £^dm = 0.27, £\ = 0.03 and Q.^ = 0.7, 
respectively. The value of the Hubble constant (in units of 100 
km/s/Mpc) is h = 0J. The initial conditions correspond to a cold 
dark matter power spectrum such that the r.m.s. matter fluctuation 
in sphere with radius of 8h~ l Mpc is G$ = 0.9. Each particle has 
a mass 6.86 x 1O 1O A _1 M0, which allows to resolve cluster-sized 
haloes by several thousand particles. The gravitational softening 
was fixed at 30h~ l kpc. 

We identified all collapsed objects from the z = output of 
the parent simulation, using a spherical overdensity criterion. From 
this catalogue we randomly selected a number of cluster-sized dark 
matter haloes and carved around each of them a spherical region 
encompassing a few virial radii. Each of these regions was then 
resampled to build new initial conditions suitable for a higher num- 
ber of particles - using on average 10 6 dark matter particles and 
the same number of gas particles. Outside this 'high resolution' 
region the number of particles was strongly reduced by interpolat- 
ing the original particles onto a coarse spherical grid with 3 de- 
gree angular resolution. This produced a number of 'low resolu- 
tion' macro particles which still give a good representation of the 
large-scale tidal field necessary to form the cluster in the high reso- 
lution region. Particle masses in the high-resolution region ranged 
from 2 x 10 9 /!~'Mq to 6 x 10 9 /! _1 M Q for dark matter, and from 

3 x IO^-'Mq to 7 x 10 8 /T 'M for gas. 

We finally evolved these new initial conditions using the pub- 
licly available code GADGET (GAlaxies with Dark matter and 
Gas intEracT; Springel et al. 2001b). This is a TREESPH code 
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Table 1. Main characteristics of the simulated clusters. Column 1: identification number. Columns 2 and 3: number of dark matter (/Ydm) and gas particles 
(iVgas) inside the virial radius R v . Columns 4 and 5: mass of the dark matter (Mdm) and gas (M gas ) components inside the virial radius. Columns 6 and 7: total 
(i.e. virial) mass M v and virial radius R v . 
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where the dark matter particles are evolved using a Barnes & 
Hut (1986) tree-code, while the collisional gas is followed us- 
ing a Smoothed Particle Hydrodynamics (SPH) approach. The 
re-simulations, which include only non-radiative hydrodynamics, 
started at redshift z; n = 35 — 50, depending on the mass resolution. 
Gravitational softening was fixed at 5h kpc (cubic spline). 

The final sample includes 17 high-resolution clusters, with 
virial masses ranging between 3.1 x 10 14 ft _1 M Q and 1.7 x 
10 A - *A/q and virial radii between 1.4 and 2.5/z~' Mpc. Virial 
radii were defined using the overdensity threshold dictated by the 
spherical top-hat model (e.g. Eke et al. 1996). 

In TableQthe main properties of the cluster sample are sum- 
marized. Our analysis is restricted to the z = outputs. A discus- 
sion of the dynamical properties of cluster substructure done using 
the results of these simulations is presented elsewhere (Tormen, 
Moscardini & Yoshida 2004). As a final comment, we stress that 
the random criterion used to select the re-simulated haloes leads 
to a sample of clusters with varying dynamical properties: at the 
present time, some are more relaxed, while others are dynamically 
perturbed. The surrounding environment can also be quite differ- 
ent: some clusters are more isolated, while others are interacting 
with the surrounding cosmic web. This means that the modelling 
we will do will be representative of an average cluster, in an aver- 
age environment and dynamical configuration. 



3 MODELLING DARK MATTER PROFILES 

As already discussed in the introduction, the underlying idea of 
this work is to build a self-consistent model for the radial profiles 
describing the dynamical properties of galaxy clusters. The general 
criterion we followed is to use at most double power-law analytic 
fits to all differential profiles. This was done in order to keep the 
number of free parameters to the minimum, and to allow easily 
integration to obtain e.g. the cluster's mass. 



3.1 Coarse-grained distribution function 

Taylor & Navarro (2001) have recently noticed that in N-body sim- 
ulations of galaxy-sized dark matter haloes the radial profile of the 
coarse-grained distribution function, f(x) = p(x) /of (x), is well fit- 
ted by a single power-law: 

7M~* a , (i) 

with a ~ —1.875. We recall that p(x) and G r (x) are the differential 
dark matter density and the radial velocity dispersion, respectively, 
at distance x = r/R v from the cluster centre, where R v is the clus- 
ter virial radius. The quantity f(x) has the dimension of a phase- 
space density, and can be thought of as the coarse-grained phase- 
space density profile of a halo. A power-law behaviour for /(x) 
is also expected in the self-similar solution for secondary infall of 
gas onto a point-mass perturber in a uniformly expanding universe 
(Bertschinger 1985). Since a single power-law is an easy profile to 
fit, we deemed it a good starting point for our investigation. 

For each cluster we computed the radial profile of f(x) versus 
the distance x from the cluster centre. We then averaged the profiles 
and fit this to a power-law relation. The result is shown in Fig.Q 
the solid line and the surrounding shaded region represent the mean 
profile with its error. 

From the figure we can notice that a single power-law fit as in 
equation Q is indeed a good description of the data at most radii. 
Our best-fitting slope, a = —1.95, is slightly steeper than that ob- 
tained by Taylor & Navarro (2001). However, at scales correspond- 
ing to logx< — 1.85, i.e. forx<0.014, the mean profile starts show- 
ing significant fluctuations from the fit. While a departure from the 
single power-law is expected for x — > (as at small radii the den- 
sity diverges, but the radial velocity dispersion goes to zero), the 
observed deviation is probably due to numerical effects. 

In order to confirm or dispute this interpretation, we com- 
puted the analogous profile for a single cluster, re-simulated with 
dark matter only, at higher and higher resolution (see Springel et 
al. 2001a). The resulting profiles are shown in Fig.Qby the three 
dotted lines: they correspond to simulations with 6 x 10 5 , 3.4 x 10 6 
and 1.9 x 10 7 particles inside the virial radius at redshift z = 0. 
These profiles, which are plotted from the radius containing 250 
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Figure 1. The coarse-grained phase-space density mean radial profile for 
dark matter in our simulations. To enhance clarity, we plot fx 2 (solid line) 
and its error (shaded region). The dashed line is our best single power-law 
fit (equation 0- The three dashed noisy lines represent the profiles of a 
single cluster re-simulated - with dark matter only - at different increasing 
resolutions (notice the different extension of the lines to the left). The ver- 
tical dashed line was taken as the limit of our model, and corresponds to 
log* = —1.85. At this distance the mean profile starts to show fluctuations 
originated by numerical artifacts. 

particles, are quite noisy because they refer to a single cluster re- 
alisation as opposed to our mean curve. However, their slope is 
quite close to that obtained from our sample, essentially confirm- 
ing our estimation of a. More importantly, the three curves show 
no systematic deviation from a single power-law at small radii: this 
seems to confirm the numerical origin of the small scale fluctua- 
tions. We therefore decided to consider all the profiles presented in 
this work only at radii larger than logx< — 1.85. Within this radius 
our clusters contain on average ~ 800 dark matter particles. This is 
a quite conservative limit, and satisfies the requirements suggested 
by recent convergence studies (Power et al. 2003). 

Notice that the normalisation of the power-law relation for 
our fit is slightly larger than that resulting from the single high- 
resolution dark matter-only simulations. This difference could ei- 
ther be due to statistical fluctuations or by the fact that our simula- 
tions contain also gas. However, the discrepancy is not really rele- 
vant for our model, as in what follows we will only use the slope oc 
of the relation. 



3.2 Velocity structure 

As a next step we will consider the dark matter velocity dispersion 
profiles. These are important not only to understand the velocity 
structure of dark matter haloes, but also to estimate the total cluster 
mass through the Jeans equation (as done in section l3~4l . 

In Fig.|2|we show the mean profiles (averaged over our clus- 
ter sample) for the radial velocity dispersion a r (x) (top panel), the 
tangential velocity dispersion O/ (x) (central panel) and the velocity 
anisotropy parameter P (x) = 1 — aj (x) /2<3 2 r (x) (bottom panel). The 
solid bands represent the error on the mean. In order to perform the 
mean over clusters of different mass, each profile was normalized 
by the quantity = GM V /R V , where R v and M v are the individ- 
ual virial radii and masses. The vertical dashed line represents the 
limit of our model, at logx = —1.85. We remind that the veloc- 



ity anisotropy (3 shows preferred directions in the velocity field. In 
particular, (3 = corresponds to an isotropic velocity field, while 
(3 < and < (3 < 1 indicate predominance of tangential and radial 
motions, respectively. 

The plots show that in general the dark matter velocity field 
in clusters is not isothermal nor isotropic. The radial and tangential 
velocity dispersions decline both at small and large radii, while the 
radial trend of (3 indicates a predominance of radial motions at large 
radii, and an almost isotropic velocity field in the inner part, a result 
consistent with previous studies (e.g. Cole & Lacey 1996; Tormen 
et al. 1997). We also find that the mean radial velocity (not shown) 
is still negative at the virial radius, i.e. infall motions are present: 
this suggests that our clusters are on average not completely re- 
laxed, as clumps of matter are still entering at z = 0. 

Since the three quantities a r , G, and |3 are not independent, 
we can derive any of them using the other two. This should also 
be true for the analytic fits; therefore we tried fitting each pair of 
profiles with double power-laws, and derived each third fit by their 
combination. The solution which best describes the simulated pro- 
files turned out to be fitting the radial velocity dispersion and the 
velocity anisotropy directly, and getting the tangential velocity dis- 
persion by combination. 

The analytical expressions of our fits (hereafter marked by a 
tilde) are the following: 

(i) radial velocity dispersion: 

ft P (x) = a r {x)/a v = o ro x°- 3 (x + x pi T 048 , (2) 

with scale radius x pi = 10~ 116 and normalization amplitude a, = 
0.67; 

(ii) velocity anisotropy (3(x): 

P(x) = p„x°- 25 (x + x / , 2 )-°- 3 , (3) 

withx K = 10 _1 1 and (3 = 1.7; 

(iii) tangential velocity dispersion: 

o f (x) = o,(x)/o v 

= d,.(x) ^2[l-p»] = C ro x 03 (x + x pi ) -°' 48 

x ^[l-pox^x + x^)- - 3 ] . (4) 

In Fig. [3] we show the mean logarithmic residuals between 
the actual profiles and our fitting formulae. Error bars represent 
the error on the mean. Notice that for radial and tangential veloc- 
ity dispersions the differences are very small and within la of the 
mean at most radii. The maximum discrepancies are found at large 
radii, but they are still below 5 per cent. The velocity anisotropy 
P(x) is equally well fitted at radii x < 0.25, while differences can 
be as large as 50 per cent at intermediate distances: 0.3<x<0.5; 
this is the range where CT,(x) is also showing the largest deviations 
from our fit. Given the non-trivial shape of the measured velocity 
anisotropy profile (bottom panel of Fig.|2j, we decided not to try 
a more complicated fit for the sake of simplicity. As we shall see 
below, the error introduced in the Jeans equation by such deviations 
is not relevant. 

3.3 Density structure 

The density profile of dark matter haloes has been widely studied 
in the last few years (e.g. Crone, Evrard & Richstone 1994; Jing et 
al. 1995; Navarro, Frenk & White 1996, 1997; Tormen, Bouchet & 
White 1997; Huss, Jain & Steinmetz 1999). There is now a general 
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Figure 2. Velocity profiles. From top to bottom, different panels refer to 
radial velocity dispersion c, (.r), tangential velocity dispersion a, (x), and 
velocity anisotropy $(x). Velocity dispersions are plotted in units of the 
virial quantity c v . The solid lines and surrounding shaded regions represent 
the mean profile obtained from our sample of simulated clusters and its 
eiTor, respectively; dashed lines represent our fitting formulae (equations l2l 
l4l see text for more details). The vertical dashed line is the limit of our 
model, set at log* = —1.85. 



Figure 3. Mean logarithmic residuals between the measured profiles and 
our fitting formulae (equations 12141 : radial velocity dispersion a r (x) (top 
panel), tangential velocity dispersion o t (x) (central panel) and velocity 
anisotropy p(jc) (bottom panel). Residuals are plotted only for radii larger 
than the fiducial limit of our fit: log,t = —1.85, as set by the analysis of 
Fig.0 
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Figure 4. Mean radial density profile for dark matter: 3?p{x) (solid line) 
with error on the mean (shaded region). The dashed and dotted lines refer 
to our fit (equation|5} and to the NFW fit (equation^, respectively. 
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Figure 5. Distribution of the concentration parameter c versus the virial 
masses M v , in units of h Mq. Each point refers to a different simulated 
galaxy cluster. 



agreement that a double power-law with a central cusp is a good 
description of simulated haloes (e.g. NFW), even though the exact 
value of the inner asymptotic slope is still debated (Moore et al. 
1998; Jing & Suto 2000; Jing 2000; Ghigna et al. 2000; Power et al. 
2003), possibly due to the presence of systematic and/or numerical 
effects. 

While it is not the purpose of this paper to enter this discus- 
sion, we would like to stress that, unlike the majority of previous 
studies, the dark matter profiles modelled in this work come from 
simulations where the effect of the hot ICM is properly taken into 
account: the presence of hot gas, even if not gravitationally domi- 
nant, can very well affect also the dark matter structure. 

A model for the dark matter density profile is straightforward 
to obtain using the fitting relations derived in the previous subsec- 
tions. If we combine the coarse-grained phase-space density profile 
/(x) (equation with the radial velocity dispersion profile 6 r (x) 
(equation |2j and solve for the density, we easily get the following 
approximate fit for the dark matter spatial density distribution: 



pW 
pb 



Po 



x(x + x pi ) 



1.5 



(5) 



with pt, the mean background density of the universe; the normali- 
sation factor po is given by 



Po 



(Wb)Av 



[(i+2x Pl )/(i+x Pl y/ 2 - 



- 2x 



1/2 



(6) 



where A v is the virial overdensity specified by the cosmological 
model; the term /j, = 0.097 is the average baryonic fraction of our 
cluster sample, measured at the virial radius; it is used to properly 
weigh the DM component in our simulated clusters. 

Notice that for simplicity we have rounded the asymptotic be- 
haviour of this fit to x -1 and x~ 2,5 at very small and very large 
scales, respectively. The exact slopes coming from the combina- 
tion of the two fitting relations for f(x) and <T r (x) would be slightly 
different, x~ 105 and x~ 2,49 respectively. However, as we will see, 
such small differences will not affect our global model. 

The mean density profile, averaged over the whole sample 
of simulated clusters, is shown in Fig. |4| As in the previous fig- 
ures, the shaded region shows the error of the mean, which is rel- 



atively small, indicating that the different clusters have quite con- 
sistent density profiles. Our model (equation [5} is represented by 
the dashed line, and agrees very well with the data, always falling 
inside the error region: only very close to the virial radius there is 
a small discrepancy. This is a consequence of the sharp change of 
slope observed in the mean profiles of both /(x) and <3 r (x), and 
cannot be improved without adding significant complexity to our 
fitting formulae. 

It is interesting to compare our relation to that originally pro- 
posed by NFW, which is largely used in the literature to model the 
halo dark matter distribution. Its expression is the following: 



Pnfw(x) 
Pb 



PO.NFW 



(x/x s )(l +x/x s ) 



(7) 



where x s is the scale parameter (inverse of the concentration pa- 
rameter: c = 1 /x s ), and corresponds to our x p , while Po.nfw is the 
normalisation factor or 'characteristic density' : 



P0,NFW : 



;i-/b)A v c 3 



(8) 



3[ln(l+c)-c/(l+c)]' 

Since pj,, A v and /j, are specified by the cosmological model, for 
any given cosmology, Po,nfw is only a function of x s , or of the 
concentration c. 

We found that the value of the NFW concentration parame- 
ter which best fits our mean density profile is cnfw ~ 6-8. This 
value is slightly lower than the mean concentration over our cluster 
sample, which is (c) = 7.13. However, the difference is well within 
the scatter in the distribution of concentrations, as shown in Fig. [5] 
where the best fitting concentration for each galaxy cluster is plot- 
ted versus the cluster virial mass. The figure also shows that there 
is no clear systematic trend with mass; consequently it is reason- 
able to fit the mean profile as representative of the full sample. The 
resulting NFW profile is shown in Fig.|4]as a dotted line. We no- 
tice that for distances 0.04<x<0.4 the NFW fit overestimates the 
mean density profiles of the simulations. 

This is more clearly visible in Fig. [6] where we plot the mean 
logarithmic residuals between the cluster profiles and the analytic 
fits (equations |5| and 0; error bars indicate the error on the mean. 
The systematic trend shown by the residuals to the NFW were al- 
ready noticed in the literature (e.g. Tormen et al. 1997), and indicate 
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Figure 6. Mean logarithmic residuals between the measured density profiles 
and the two fitting relations. Top and bottom panels refer to equations 151 
and 0, respectively. Error bars represent the error on the mean. 



that the NFW fit overestimates the actual profiles in the simulations 
by up to 15 per cent at x ~ 0. 1, with a rms residual of the mean of 
12.1 per cent. This behaviour is strongly reduced if the relation pro- 
posed in this paper is used: in this case residuals are always within 
the error bars, except at distances close to the virial radius, and 
have a rms of 6.0 per cent. This result is not necessarily at odd with 
the literature; we wish to remind that the NFW fit was proposed 
to model the density structure of isolated and relaxed dark matter 
haloes in dark-matter-only simulations. The fact that our clusters 
are not all isolated and relaxed, and the presence of a hot gas com- 
ponent, can both play a role in modifying the dark matter profile. 
We will come back to this issue in the discussion section. 



3.4 Mass Estimates 

Given the previous relations, the cluster mass inside any given ra- 
dius can be obtained either by direct integration of the density pro- 
file, or by applying the Jeans equation to the system. 

Integration of the density profile (equation |5} gives immedi- 
ately 



Mix) 



M(<x) 

My 



6po 



x + 2x D 



-x p y/2 



■2x 



1/2 



(9) 



The corresponding gravitational potential profile has also an ana- 
lytic expression: 



= -16jiGp 



(x+ Xp y/ 2 -x 1 / 2 



(10) 



For comparison, we also report the corresponding relations 
obtained by integrating the NFW profile (equation^: 



Mnfw(<x) _,. in{l+cx)-cx/(l+cx) 
~M y ~ ( ~ M ln(l+c)-c/(l+c) 

for the mass, and 

47tGpo,NFW , , , , 

*NFW (x) = 3 In ( 1 + cx) , 



(11) 



(12) 



for the gravitational potential. 

However, it is more interesting and useful to estimate the mass 



profile through the Jeans equation which, for a spherical static sys- 
tem, can be expressed as the ratio between the estimated and true 
mass: 

xR v a 2 r {x) 



M h (<x) = — 



dlnp(x) dlna r (x) 



rflnx 



d\nx 



+ 2PW 



(13) 



Tormen et al. (1997) and Thomas et al. (1998) have shown that 
cluster-sized dark matter haloes are indeed reasonably modelled by 
this equation. Observed galaxy clusters have also been described 
by such an equation (e.g. the CNOC cluster survey, Carlberg et al. 
1997) with results in agreement with those coming from simula- 
tions. We stress that we are not imposing the Jeans equation as a 
constraint to our model: we are using it to test the self-consistency 
of our density and velocity fitting formulae, assuming that the 
galaxy clusters obey it. Using the Jeans estimate allows in principle 
to also weigh the different contributions to the total estimate, and 
to measure the errors made when one or more terms in the equation 
are ignored for lack of information on the system. 

The performance of the different mass estimators is shown in 
Fig.0 where we plot the ratio between the actual mass profile (i.e. 
the mean mass profile of our simulated clusters) and the mass pro- 
file obtained by either integrating the density profile or by inserting 
the analytic fits into the Jeans equation 1 1 3i . The profiles from sim- 
ulations were multiplied by a factor (1 +e 2 /* 2 ) 3 ' 2 , where e is the 
Plummer equivalent of the actual spline softening, to deconvolve 
the effect of gravitational softening and thus provide a more ac- 
curate comparison at small radii. Curves labeled 1 and 2 refer to 
estimates based on the integrals of the density models, M{x) (equa- 
tion|5J and Mnfw(< x) (equation lllk respectively; curve 3 refers 
to the mass obtained from the Jeans equation 1 1 31 using the density 
and velocity profiles given in equations J2j, {3J and |5). 

Notice a very good agreement between the actual mass and 
that derived by integrating our density profile (equation [5j- The 
maximum discrepancy is found at small radii (log* ~ —1.85) and 
is of the order of 10 per cent. The analogous estimate from the 
NFW density profiles shows larger discrepancies: roughly 40 per 
cent at logx ~ —1.85 and 10 per cent at log* ~ — 1. On the other 
hand, both models accurately reproduce the mass profile close to 
the virial radius, where the errors are of the order of 1-2 per cent. 

The mass estimate obtained by the Jeans equation (curve 3) 
is also very good; its error is always smaller than 10 per cent, sug- 
gesting that the model we are proposing is indeed dynamically self- 
consistent. The main discrepancy is found at distances close to the 
virial radius; we have checked that this error originates mainly from 
the incorrect slope of our density profile close to R v . 



4 MODELLING GAS PROFILES 
4.1 Density structure 

In analogy to what has been done in the previous section for the 
dark matter profiles, in this section we will apply a similar analysis 
to the radial profiles of the hot gas component. Let us start from the 
density. Fig. [S] compares the profiles of gas and dark matter, nor- 
malised so that they match at the virial radius. The plot shows that 
the profiles are similar roughly for x>0.06, while the gas profile 
becomes flatter in the internal region. The different internal slope 
of the two profiles agrees with previous analyses, coming from sim- 
ulations at smaller resolution (e.g. Navarro, Frenk & White 1995; 
Eke, Navarro & Frenk 1998; Frenk et al. 1999; Pearce et al. 2000; 
Lewis et al. 2000; Yoshikawa, Jing & Suto 2000). Recent analyti- 
cal models of the ICM in galaxy clusters (Komatsu & Seljak 2001) 
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Figure 7. Ratio between the true mass M Uuc (< x) and the one derived from 
different estimators. Different curves refer to M(x) (equation|9] solid curve 
1), Mnfw(< x) (equation ll II dashed curve 2); and M E {< x) (equation ll3l 
dotted line 3). 
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Figure 8. Dark matter (dashed curve) and gas (solid curve) density radial 
profiles. The two curves are normalized to have the same value at the virial 
radius. The vertical line represents the limit we assumed for our model, at 
log.v = -1.85. 



also propose a gas density profile with a flatter inner part. In our 
case, however, the similarity extends to inner radii than previously 
found. 

The similarity of the DM and gas profiles in the outer part of 
the cluster suggested us a first naive approach: that of fitting the 
gas density profile with the same functional form used for the dark 
matter, with a different concentration value, assuming of course a 
different normalization based on the average baryon fraction. The 
best-fitting NFW profile for the mean gas density turned out to have 
a concentration c = 4.3. However, the discrepancy (not shown) be- 
tween this model and the mean profile in the simulations resulted 
quite large, the fit having too steep a slope in the inner part, where 
the gas profiles becomes flatter than x . 

We therefore relaxed the constraint on the inner slope, and 
tried to fit the gas density profile with a double power-law having 
the same outer slope as the DM profile, but a different asymptotic 



behaviour at small radii. The best fitting result turned out to be the 
following analytic expression: 



PW 
Pb 



Po 



where x„ = 0.04; the normalisation po is given by 



Po ; 



3 



16^/3 + 40^/3 + 10x p + 2 164' 



1/2 



.1)2.5 



(14) 



(15) 



This fit tends to a constant value in the limit x — > and de- 
creases as p(x) «: x -2 ' 5 at large radii. The upper panel of Fig. [5] 
compares the proposed fit (dashed curve) to the average profile for 
the gas density obtained from our simulated clusters (solid curve): 
the agreement is always very good and the fitting relation falls well 
inside the shaded region, which displays the l-o" error of the mean. 
Notice that the results for our simulated clusters are shown starting 
from the distance containing at least 250 gas particles, while the 
fitting relations are computing for logx > —1.85, only. 

A possible different approach is based on the so-called (3- 
model (Cavaliere & Fusco-Femiano 1976, 1978) which is the den- 
sity profile most commonly used to describe the gas distribution 
in observed galaxy clusters. Here the original underlying assump- 
tion is that the gas density is proportional to the galaxy density to a 
given power (3. Since the galaxy distribution was assumed to follow 
a King profile, the resulting gas density profile is given by 



PpW 



Pop 



(,v2 



\W2 



(16) 



where x c is the core radius (in units of the virial radius) and the 
external slope is given by 3(3. However, more recent analyses have 
shown that the galaxy distribution traces the total mass distribu- 
tion, which is well described by a NFW profile (e.g. Carlberg et al. 
1997). Moreover, equation 1 16i cannot fit - at the same time - the 
central and external parts of the cluster profiles, as traced by SZ and 
X-ray observations. Finally, previous numerical works have already 
shown that the |3-model is unfit to describe the results of (lower res- 
olution) hydrodynamical simulations (e.g. Navarro, Frenk & White 
1995; Bartelmann & Steinmetz 1996). Still, it is useful to test the 
(3-profile against the high-resolution simulations presented in this 
work. We computed the two parameters x c and (3 by best fitting 
the mean profile of the gas density, and obtained x c = 10~ 1,53 and 
(3 = 0.7. The corresponding (3-profile is shown as a dotted curve in 
the upper panel of Fig.|9j It is clear that there is not agreement with 
the mean result of our simulations, being too flat at large radii. In 
order to quantify the goodness of fit for the two models, in the lower 
panel of the same figure we show the logarithmic residuals between 
the gas density of simulated clusters and our fit (top) or the (3-model 
(bottom). While in the first case the ratio is very close to unity for 
x<0.8, the B-model has a 10 per cent deviation for 0.1<x<0.3 
and overestimates the density by 25 per cent at the virial radius. 



4.2 Velocity structure 

We next consider the gas velocity dispersion profiles, normalized 
by the virial velocity dispersion o v as we did for the dark mat- 
ter. Results are shown in Fig. 1101 where the solid curve and the 
shaded region represent, as usual, the mean values from the simula- 
tions, and the error on the mean. Only for the profile of the velocity 
anisotropy (3(x) we make use of the median instead of the average: 
since (3(x) usually fluctuates around zero, an average curve would 



A dynamical model for the distribution of dark matter and gas in galaxy clusters 9 



10 




0.01 



0.1 

x = r/R 



1 = 



0.9 



0.7 



3 °- 9 r 

cT 0.8 | 



0.7 



H I I I I I I : 



t iWiii - 



0.01 



0.1 

= r/R v 



Figure 9. Gas density profile. Top panel: average radial profile px 2 (solid 
curve) and its error (shaded region). The dashed and dotted curves represent 
our model (equation 1141 and a P-model (equation 1161 , respectively. The 
vertical line indicates the limit of our fit, taken at log* = —1.85. Bottom 
panel: logarithmic residuals between the simulated gas density profiles and 
the analytic models. The upper and lower parts of the plot refer to our model 
and to the p-model, respectively. The points and error bars represent the 
mean ratio and their errors. 



be strongly affected by numerical fluctuations, and would bias the 
mean curve away from zero. 

We observed, as expected, that the gas velocity in Fig. llOl be- 
haves very differently from the DM of Fig. [2] due to the different 
kind of forces acting on the two components. Both the radial and 
tangential gas velocity dispersions are increasing functions of the 
distance x, even if a r (x) presents a central flattening. The veloc- 
ity anisotropy (3(x) is marginally negative in the intermediate part 
of the clusters (x ~ 0.1), denoting a slight predominance of tan- 
gential motions. Finally, the gas radial velocity (not shown) is not 
completely zero, but slightly negative and larger at large radii, sug- 
gesting a residual net infall. 

Our best analytic fits to the curves for a r and a t are also shown 
in Fig. 1101 indicated by dashed lines; their expressions are: 



6>(x) = o>(x)/o\. = o- ro 
where o", = 0.45 and x„ 



(17) 
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Figure 10. Radial velocity profiles for the gas: radial velocity dispersion 
<S r (x) and tangential velocity dispersion o,(,t) (top and central panels, re- 
spectively; the curves are normalized as in Fig.|2j and velocity anisotropy 
P(x) (bottom panel). In the upper and central panels solid lines and shaded 
regions are the mean values obtained from our sample of simulated clusters 
and their errors, respectively; only for p we prefer to show the median val- 
ues and its errors. Dashed lines show our fitting formulae (equations ll7H9l . 
The vertical dashed line represents the limit of our model, log* = —1.85. 
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Figure 11. Logarithmic residuals between the simulated gas velocity pro- 
files and our models (equations 1171181 . Radial velocity dispersion O r (x) 
and tangential velocity dispersion a, (.v) are shown in the upper and bottom 
panels, respectively. The points and error bars represent the mean ratios and 
the errors on the mean. The vertical dashed line represents the limit of our 
model, logx= —1.85. 



C t (x) =0/(x)/0 v = 



.!).<-, 



(x P2 +x)° 



(18) 



where a to = 0.54 and x pi = 10 



The logarithmic residuals between the mean simulated profiles 
and our models are plotted in Fig. ll II with mean values and error 
indicated as usual. The error bars are quite large, showing a large 
variance in our cluster sample, but the results are always compatible 
with unity, denoting a good agreement between simulated clusters 
and models. 

In order to write a model for the velocity anisotropy p(x) we 
decided to simply combine the two previous fits for the velocity 
dispersions a t (x) and a r (x): 



P(x) = l- 



26 2 (a 



1 - ■ 



a 2 x 1 - 2 



2a r ( x p, +x){Xp i +x) 



(19) 



The resulting curve, displayed as dashed line in the correspond- 
ing panel of Fig. llQI well reproduces the behaviour on the median 
value, showing the self-consistency of our relations. 



4.3 Temperature structure 

The radial temperature profiles of galaxy clusters, as estimated 
from X-ray observations, still show large uncertainties. Only re- 
cently, thanks to spectroscopic data coming from the Beppo-SAX 
and Chandra satellites (see, e.g., Allen, Schmidt & Fabian 2001; 
De Grandi & Molendi 2002; Ettori, De Grandi & Molendi 2002a; 
Ettori et al. 2002b; Johnstone et al. 2002), the complexity of the 
thermal structure inside galaxy clusters has been revealed. In par- 
ticular, observations suggest the presence of a central isothermal 
region followed by a smooth decline towards the centre. In order to 
reproduce this feature, different recipes based on the inclusion of 
radiative cooling and feedback heating have been unsuccessfully 
attempted in numerical simulations (Lewis et al. 2000; Loken et 
al. 2002; Muanwong et al. 2002; Valdarnini 2003; Tornatore et al. 
2003; Borgani et al. 2004). We remind that these processes are not 
considered in our non-radiative hydrodynamical simulations; how- 
ever, they influence only the very central part of galaxy clusters, 
and so in practice they do not affect the model we propose here. 
The behaviour of the temperature profile is also observationally not 
well determined at large radii and shows large variation from ob- 
ject to object. However, there are strong indications that the value 
of the temperature at the virial radius is approximately a factor of 2 
smaller than at the centre. This fact in turn causes errors in the mass 
estimates obtained using the equation of the hydrostatic equilib- 
rium. Lacking a reliable measure of the radial temperature gradient 
in clusters based on X-ray observations, it is highly useful to have 
some indications from numerical simulations. For these reasons we 
devoted great care in modelling the mean temperature profile. 

In the upper panel of Fig. 1121 we show the average tempera- 
ture profile normalized by the estimate of the virial temperature, 
T v = (G/.im p M v )/(k n R v ). In this relation G is the gravity constant, 
jj = 0.59 is the mean molecular weight, m p is the proton mass and 
k\, is the Boltzmann constant. The figure clearly shows a nearly 
isothermal inner region, out to x ~ 0.2, with only a hint of positive 
gradient. At larger radii the temperature starts to decrease, reaching 
at the virial radius a value about 60 per cent of the central one. The 
profile we find is in good agreement with previous numerical works 
(Navarro et al. 1995; Eke et al. 1998; Bryan & Norman 1998; Frenk 
et al. 1999; Thomas et al. 2001; see, however, Loken et al. 2002; 
AYMG03) and also with the temperature profile found from the 
analysis of the Beppo-SAX data by De Grandi & Molendi (2002). 
However, we remark that the much higher numerical resolution of 
these simulations makes us more confident in the result and allows 
us to model the simulated profile with more accuracy. 

The sudden temperature drop at x ~ 0.2 cannot be fit it by a 
double power-law as those used so far. We model this sharp transi- 
tion using instead the following expression: 



f(x) = T(x)/T v 



T x« 



(x4 +4)0.13' 



(20) 



where Tq = 0.255, x p = 10 U M ; the internal and external logarith 



p> 
10 -0.51. 

mic slopes are approximately and —0.5, respectively. The bottom 
panel of Fig. ll2l shows the logarithmic residuals between the simu- 
lated profiles and the previous relation. The mean ratio is very close 
to unity, with no systematic trend, confirming the accuracy of our 
fit. 



4.4 Mass Estimates 

In analogy to the dark matter, we can derive a first estimate of the 
mass of hot gas inside any given radius by integrating the gas den- 
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Figure 12. Gas temperature profile. Top panel: mean radial profile T(x) 
(normalised to the virial estimate T v ; solid line) and its error (shaded re- 
gion). The dashed line refers to our model (equation l201 . The vertical line 
represents the limit of our model, logx = — 1 .85. Bottom panel: logarithmic 
residuals between the simulated cluster temperature profile and the analytic 
model. The points and error bars represent the mean ratio and the errors on 
the mean, respectively. 



sity profile (equation ll4t . We obtain 



A*icmM 



Micm(<*) 



My 



lxp/3- 



6pp 



20xx 2 p /3 + 5x 2 x p - 



8x 



1/2 



12.5 



(21) 



The short-dashed curve labeled 6 in Fig. II 31 is the ratio between 
the actual gas mass and the mass obtained from equation <2U : the 
agreement between actual and modelled mass profiles is excellent, 
with differences smaller than 5 per cent at all radii x > 0.02. The 
expression in equation <2 1 1 will be used below to discuss the bary- 
onic fraction of our clusters. 

We now have all the ingredients to discuss other mass esti- 
mates using the gas profiles. One of the standard assumption used 
to derive cluster masses from X-ray data is that the system is in hy- 
drostatic equilibrium: the results of our high-resolution simulations 
can give us some indications on the accuracy of this assumption, 
allowing us to evaluate the size of the typical errors made in clus- 



ter mass estimation when complete dynamical information are not 
available. 

Assuming a spherical and static gravitational potential, the 
usual equation for hydrostatic equilibrium can be formally obtained 
from the Jeans equation (equation ll3> by equating the gas internal 
energy to the dark matter velocity dispersion (T r DM> 



J /\DM 



(22) 



lim p 

and by assuming an isotropic velocity field: P(x) = 0. One then 
obtains the following mass estimate 



M E {<x) = - 



xR v k^T(x) 
Gfmip 



dlnp(x) d\nT{x) 



dlnx 



dlnx 



(23) 



However, the hypotheses underlying equation I23i are too restric- 
tive for at least two reasons. First, our simulations show that the 
mean velocity anisotropy P(x) is not zero inside the virial radius 
(see Fig. UOt . Second, the replacement of the velocity dispersion 
with the temperature is not fully justified, due to the residual bulk 
motions shown in the same figure. 

In order to have a more complete description of the gas be- 
haviour, it is more convenient to start directly from the force equa- 
tion, and to keep all terms coming both from isotropic pressure and 
from anisotropic velocity dispersion. For a system in a spherical 
and symmetric potential <fr(x) = GM(< x)/x, this gas dynamical 
equilibrium equation can be written as: 



d&{x) _ 1 dP{x) 1 d[p(x)a 2 (x)] 



dx 



p(x) dx p(x) 



dx 



-2p(x 



0?W 



(24) 



By substituting the expression for the pressure, P(x) = 
p(x)h D T(x)/iAm p , it is possible to obtain a new mass estimator 
M E (< x) = M(< x)/M y for the total mass inside x traced by the 
gas distribution: 



xR v k\,T(x) 

GfMp 
X i? v O 2 (x) 



dlnp(x) dlnr(x) 
dlnx dlnx 



dlnp(x) d\n<3 2 (x) 



dlnx 



dlnx 



■ (25) 



We then used this complete equation, together with some simplifi- 
cations of it, to estimate the cluster mass inside x using the radial 
profiles we proposed for the various quantities. In particular we 
considered the following possibilities: 



(i) a P-model isothermal sphere, giving: 
xR v k h T(<x) dlnpp(x) 



Mf (<x) = — 



dlnx 



(26) 



where T(< x) is the actual mean temperature inside x and pp(x) is 
the P-model defined by equation i 161 : 

(ii) an isothermal sphere with density profile given by our best 
fit to the simulations, equation i 14t . giving: 



<x =- 



xR v k^T(< x) d In p (x) 



Gfjm p 



dlnx 



(27) 



(iii) the usual hydrostatic equilibrium equation <23t with a P- 
model gas density profile (equation ll6t and our best fitting relation 
for the temperature profile (equation l20> : 



Mf(<x 



-XMyf (X) 



dlnpp(x) dlnf (x) 



(28) 



dlnx dlnx 
(iv) the same equation of hydrostatic equilibrium I23i . where 
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Figure 13. Ratio between the true mass M tlll e and the mass derived from 
different estimators M est presented in the text. The curves refer to Mf (equa- 
tion^ dotted line), Mf (equation^] long-dashed line), Mf (equationl29l 
short-dashed line), Mf (equation l28l dotted-short-dashed line), Mf (equa- 
tion l30l solid line). Curve labeled 6 represents the ratio between the actual 
gas mass and the mass obtained from equation <2 1 1 : dotted-long-dashed 
line. 



the density and temperature profiles are our best fit equations 1 141 
and |20}: 



Mf(<x) 



-xM v f(x) 



dlnp(x) d\nf(x) 



dlnx dlnx 
(v) the more general model described by equation 1251 : 
dlnp(x) dlnf(x) 



(29) 



Ml(<x) 



—xM v T(x) 



+xM v 5 2 r (x) 



dlnx dlnx 
dlnp(x) dlna^(x 
dlnx dlnx 



+ 2p(x 



,(30) 



where all the profiles are given by our best fits to the simulations, 
equations (HJ, H9. (13 and (191 ; 

As done for the dark matter component in Fig. we dis- 
cuss the performance of the different mass estimators by plotting in 
Fig. ll3l the ratio between the true mass profile, i.e. the average mass 
profile of our simulated clusters, and the previous analytic relations. 
Curves labeled 1 and 2 refer to estimates based on the isothermal 
sphere model (Mf using equation [26] and Mj using equation |2*7V 
curves labeled 3 and 4 are for a standard hydrostatic equilibrium 
model using for the temperature the new fit we propose and for the 
gas density a P-model (M^, equation |28K or using for temperature 
and density our new fits (Mf , equation l29t ; curve labeled 5 is for 
the complete 'hydrodynamical equilibrium' model which includes 
the support from residual gas motions, and uses all the analytical 
relations we derived (M^, equation l30l . 

From the figure we notice that the p-models (curves 1 and 3) 
underestimate the true cluster mass at most radii, the disagreement 
being worst at large radii, where the isothermal (3-model (curve 1) 
is off by up to 30 per cent, and the non-isothermal one (curve 3) is 
off by up to 50 per cent. This is in agreement with previous results 
in the literature (see e.g. Muanwong et al. 2002). Using our density 
profile for the gas (curve 2 and 4) improves the estimate at all radii 
of interest (x > 0.2), reducing the mass error at the virial radius 
to 15 per cent (isothermal model, curve 2) and to roughly 30 per 
cent (non-isothermal model, curve 4). The complete model (curve 



5) is the most accurate at all radii, with an error of 5 to 10 per cent 
at x > 0.2. Comparison of curves 4 and 5 shows that the pressure 
contribution from bulk gas motions is significant at all radii x > 0. 1 . 

It is also interesting to see that the two isothermal estimators 
(curves 1 and 2) on average fare better than the non-isothermal ones 
(curves 3 and 4) at all x > 0.2. The reason for this counterintuitive 
result lies in the different temperature used in the isothermal mod- 
els (average temperature within each radius) compared to the non- 
isothermal ones (local temperature at each radius). By chance this 
difference cancels part of the error made by ignoring the gas mo- 
tion in the outer cluster regions, and effectively makes isothermal 
estimators more accurate than non-isothermal ones. This is an inter- 
esting result: if it were confirmed also when projection effects and 
observational uncertainties are considered, it could lead to more 
accurate mass estimations from X-ray observations. 

We also note that the incomplete estimators (curves 1 to 
4) generally underestimate the actual cluster mass, but a simple 
isothermal model with the appropriate gas density profile (curve 
2) is still accurate to the 15 per cent level. 



4.5 Baryonic Fraction 

A straightforward consequence of our model is an analytic model 
for the baryon fraction radial profile, obtained by using the expres- 
sions for the dark matter and gas masses (eq uations l9l and l2Tl : 



/bar(<-*) 



M„ 



<X 



M g - ds (<x)+M DM (<x) 



(31) 



Notice that / bar (x < 1) = f b = 0.097. 

The previous expression is compared to the actual baryon 
fraction in Fig. 1141 The analytic profiles predict a relation (solid 
curve) that is practically indistinguishable from the actual profile 
for x > 0.1, while at smaller radii deviations grows, reaching a 15 
per cent overestimate at x = 0.02. However, this can still be consid- 
ered an excellent result of our analytic fits. 



5 DISCUSSION 

5.1 Effect of ICM on the density profile 

The results of Section l3~3l on dark matter density profiles indicate 
that the NFW fit has systematic residuals at intermediate radii with 
respect to the actual mean dark matter profile of our simulated 
haloes. Although this was already shown by Tormen et al. (1997) 
on a set of dark matter-only simulations, it is natural to ask more 
generally whether and how much the presence of a hot gas com- 
ponent affects the dark matter distribution. All the SPH clusters 
studied in this work were also simulated with dark matter only, us- 
ing otherwise the same initial conditions and identical parameters 
for the evolutionary code: we can thus directly compare the density 
profiles in the two instances. 

We would like to stress that the NFW profile was originally in- 
troduced as a fit for the dark matter profile of isolated and relaxed 
dark matter haloes in pure N-body simulations, and later it has been 
used to fit the total profiles of observed clusters. In both cases the 
fitted profiles represent the total density distribution, not only the 
dark matter one. Therefore, we think it is here more appropriate to 
compare the dark matter-only profiles of pure N-body clusters to 
the total (i.e. dark matter plus hot gas) profiles of our SPH clusters. 
These are shown in Fig. 1151 the short- and long-dashed curves rep- 
resent the mean dark matter profiles in the SPH runs and in the pure 
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Figure 14. Baryonic fraction. Top panel: radial profile of the baryonic frac- 
tion. Solid curve and shaded region indicate the mean and its error for the 
cluster sample; the dashed curve is the baryonic fraction predicted by in- 
tegrating the analytic profiles proposed in this paper (equation l3U . Bottom 
panel: logarithmic residuals between the actual and predicted baryonic frac- 
tions. 



N-body companion runs, respectively, while the dotted curve is the 
gas density (offset by a factor f^ 1 for clarity), and the solid curve 
is the sum of dark matter and gas profiles in the SPH runs. Compar- 
ison of the two total profiles (long-dashed and solid curves) shows 
that the effect of gas is to slightly concentrate the total distribution. 
An NFW fit on the long-dashed curve (total N-body density) gives 
a mean concentration cj^ ~ 6.0, while a fit on the solid curve (to- 
tal SPH density) gives a concentration csph ~ 6-5. Thus the effect 
of adding a hot gas component shifts concentrations up by roughly 
10 per cent, an amount comparable to the baryonic fraction of our 
SPH runs. 

The explanation of this result might be found in the different 
energy content of DM and gas. As first pointed out by Navarro, 
Frenk & White (1995), during the accretion of matter onto a clus- 
ter, the collisional nature of the ICM delays the gas infall compared 
to the dark matter. Due to this spatial lag, the gas feels a slightly 
stronger gravitational field and so acquires extra energy at the ex- 
penses of the dark component. This fact was explicitly demon- 



strated by Pearce, Thomas & Couchman (1994) and by Navarro 
et al. (1995), but can be also observed by comparing the energy 
content of the two species through the ratio 



P. 



dyn 



k h T/pim p + al 



(32) 



where (7dm an d °gas are the one-dimensional velocity dispersions 
for the dark matter and gas, respectively, and T is the gas tem- 
perature. This parameter differs from the usual (3 sp ec for the extra 
term <3g as , which accounts for the gas bulk energy. In Fie. 1 161 we 
show the radial profiles of the average value of pd yn m spheres of 
given radius, Pdyn ( < *)■ Solid curve and band indicate mean over 
the cluster sample and its error. This quantity illustrates the en- 
ergy ratio of dark matter and gas inside a sphere of any radius x. It 
clearly shows that, at radii larger than about x = 0.1, there is energy 
equipartition between the two species. Instead, at smaller radii, the 
gas has a larger energy content than its collisionless companion. 
The same radius is also close to where the similarity of the gas and 
dark matter density profiles breaks down, as shown in Fig. [8] the 
extra energy acquired by the gas makes it further expand, resulting 
in a flatter profile than that of the DM. For reference, the curve with 
hatched band is the profile of p spe c(< x), and shows that neglecting 
the gas bulk energy hides the existing energy equipartition at large 
radii. 

The result of Fig. ll5l seems to indicate that the energy lost by 
the dark matter allows the halo to further compress and concentrate 
compared to a pure collisionless collapse. We may speculate that, 
even if not large per se, this effect could lead to observable differ- 
ences in the properties of galaxy clusters at small scales, e.g. with 
respect to gravitational lensing effects. 

How well are these profiles fitted by the NFW formula? The 
rms residual of the mean profile for the pure N-body simulations is 
16.7 per cent, comparable to previous findings (e.g. Tormen et al. 
1997); the residuals of the mean DM profile in SPH runs is 12.1 per 
cent; the residual on the total (DM plus ICM) mean profile in the 
SPH runs is 13.3 per cent: therefore, besides changing the concen- 
tration, adding the hot gas component slightly improves the NFW 
fit compared to the pure N-body runs. This seems to indicate that 
the departures of the DM density from an NFW model are not due 
to the addition of gas (which in fact slightly improves the fit), but 
rather to the fact that our clusters were chosen with no specific cri- 
terion on their dynamical status: they are not particularly relaxed, 
nor necessarily isolated. Just for reference, if we try to fit the same 
profiles with our relation for the dark matter (equation |5J we ob- 
tain the following rms residuals of the mean curve: 9.4 per cent, 
6.0 per cent and 4.6 per cent for DM in N-body runs, DM in SPH 
runs and total density in SPH runs, respectively. Therefore, the fit 
(equation|5J constitutes an improvement over the NFW in all cases. 



5.2 Dependence on the dynamical and environmental status 

As already said, we selected our sample in a random way, with no 
constrain on the dynamics or environment of the objects. On the 
other hand, the NFW model is known to apply better to clusters 
which are relaxed and virialized. It is therefore important to inves- 
tigate how the fitting relations we propose in this paper depend on 
the dynamical and environmental status. 

To this aim, we first created a subsample excluding the sys- 
tems which show at z = evidences of a major merging event. Fol- 
lowing Tormen et al. (1997), we defined as perturbed those clus- 
ters for which the total mass, the dark matter velocity dispersion 
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Figure 15. Comparison of density profiles in different simulations. The 
long-dashed curve is the average dark matter profile in the pure N-body 
runs; the short-dashed curve is the average dark matter profile in SPH sim- 
ulations; the dotted curve is the average gas profile in the SPH simulations, 
moved up by a factor fC ' , for clarity; the solid curve is the total (dark matter 
plus gas) profile in SPH simulations. 





Figure 17. Dark matter radial density profile for the subsample of 5 isolated 
clusters. The average and the corresponding error are shown by the solid 
line and shaded region. Dashed and dotted lines refer to our fit and to the 
NFW fit, respectively. 



particular, this is true for the dark matter density profile: this lets 
us conclude that the difference between the NFW and our density 
profiles seems not to be due to the presence of merging systems in 
the original sample. 

A slightly different approach to select well-behaved objects 
is considering the environment of each system: following Lanzoni 
(2000), we defined as isolated those clusters (with virial mass My 
and radius Ry) for which no companion with mass larger than 1/100 
My was found inside a sphere of radius 2 Ry. In this way we se- 
lected, from the original sample of 17, a subset of only 5 isolated 
clusters (which is also a subset of the 12 relaxed ones). The re- 
sulting average profiles are similar to the previous ones, but now 
the DM density profile (Fig. 1171 shows a better agreement with the 
NFW fit: in particular, the outer slope is somewhat steeper, and the 
profile at intermediate region is more curved, resulting in a reduc- 
tion of the residuals. 

This fact suggests that the average density DM profile for our 
sample deviates from the NFW fit mainly because of the inclusion 
of non-isolated objects. 



Figure 16. Cumulative profile of pd yll . The curve and solid band show the 
mean value (and its error) for equation 1321 calculated in spheres of radius 
x. The curve with hatched band refers to p s „ ec , i.e. it is obtained from the 
same equation, but neglecting the gas bulk energy. 



and the gas temperature show a sudden significant increase (say at 
least 20 percent in mass) between the two last available outputs. 
We found 5 such objects: excluding them we obtain a subsample of 
12 relaxed clusters, for which we computed the average profiles of 
all relevant quantities. We found that the resulting average profiles 
are statistically identical to those coming from the full sample; the 
only exception is the radial velocity, which is much closer to zero at 
the virial radius for the relaxed sample. This is consistent with the 
fact that selecting relaxed objects we tend to exclude systems with 
significant infalling of matter. More importantly for our considera- 
tions, we found that all the fitting relations described in the previous 
sections give again a proper description also for this subsample. In 



5.3 Comparison with Ascasibar et al. (2003) 

The recent paper by AYMG03 has presented an analysis similar to 
ours, applied to a sample of 15 simulated groups and poor clusters. 
Their approach is different in spirit, as they assume that the gas is in 
hydrostatic equilibrium and that it follows a polytropic equation of 
state, whereas we verify our model with the hydrostatic equilibrium 
equation, and do not assume any equation of state. Nevertheless, 
we can compare some of their results with ours, once the different 
definitions of cluster radii are taken into account. 

In agreement with AYMG03 we also find that the p-model 
is not a good fit to the gas profile of simulated clusters. Instead, 
we do not agree on the temperature profiles. AYMG03 claim that 
a central isothermal structure is a numerical artifact of standard 
SPH schemes which do not implement entropy-conservation, and 
in fact we used the public version of GADGET, which does not use 
the entropy-conservation algorithm. However, AYMG03 claim that 
differences between the two SPH schemes are significant only at 
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scales larger than ~ 20 kpc for their Cluster A, and indeed their 
Fig. 1 indicates that above ~ 0.03i? v the gas density and tempera- 
ture profiles in the two SPH implementations are statistically indis- 
tinguishable. On the other hand, our simulations exhibit isothermal 
cores out to x ~ 0.2, a much larger scale. Although we do not have a 
clear explanation for this discrepancy, it could be that the disagree- 
ment is at least enhanced by the different sample used (a sample 
of rich clusters versus a sample of groups and poor clusters) and 
by the different dynamical status (if the sample of AYMG03 is less 
dynamically relaxed than ours, then the central temperature of their 
systems would be on average perturbed and biased towards higher 
temperatures). We are currently creating and analysing a sample 
of clusters at even higher resolution, simulated with the entropy- 
conserving SPH scheme; the analysis of this new sample should 
provide a definitive answer to this question. 

As a consequence of having an isothermal core, we find that 
a single polytrope is not a good description of the gas in our clus- 
ters. While at radii x > 0.2 the polytropic index y is close to 1 .2 
(the value quoted by AYMG03), at smaller radii y drops to unity. 
AYMG03 plot their polytropic index only down to 0.2/?20Ch roughly 
corresponding to x = 0.15, so we cannot make a direct comparison 
at small radii. 



6 CONCLUSIONS 

We have proposed a set of analytical fits which accurately describe 
the mean radial properties of dark matter and diffuse gas in simu- 
lated galaxy clusters. Our relations have been calibrated on a sam- 
ple of 17 high-resolution SPH clusters, chosen at random from a 
large parent cosmological simulation. The model provides a self- 
consistent description of these systems, and can be used to write 
accurate mass estimates by mean of the Jeans equation for the 
dark matter distribution and of a similar hydrodynamical equation 
for the diffuse gas distribution. For dark matter we propose a new 
model for the density profile, obtained combining phase-space in- 
formation and radial velocity dispersion; the resulting profile de- 
scribes our simulations more accurately than the NFW model at 
r < 0.7R V . This model differs from the NFW one in that the asymp- 
totic density slope at large radii is —2.5 instead of —3. 

More importantly, we introduce simple new fitting formulae 
for the DM radial and tangential velocity dispersion profiles and 
for the velocity anisotropy profile. Radial and tangential velocity 
dispersions profiles have a maximum value at intermediate radii 
and decrease both towards the cluster centre and outskirts. Velocity 
anisotropy is an increasing function of radius: orbits are isotropic 
in the centre and become increasingly radial orbits at larger radii. 

These profiles are tested against the Jeans equation, and are 
shown to give a self-consistent dynamical model, useful also for 
mass estimates. We find that the Jeans equation gives a realistic 
estimate of the mass (with error typically less than 10 per cent). 
We stress that the focus of our analysis of the dark matter profile 
is not to propose a replacement for the NFW fit, but to extend the 
model to all the dark matter profiles relevant for studying the cluster 
internal structure and dynamics. 

We performed a similar study for the gas profiles, introducing 
new analytic fits for gas density, temperature and velocity disper- 
sions. The average gas density profile of our sample is well de- 
scribed by a double power-law, similar to the dark matter density at 
r > 0.06/? v , and approaching a constant density towards the centre. 
The average temperature profile shows an isothermal core extend- 



ing to r ~ 0.2R V , followed by a steep decrease that reaches a factor 
two lower around the virial radius. 

The gas velocity is not zero inside the clusters: residual mo- 
tions become smaller towards the centre - showing that the infalling 
gas slows down as it sinks towards the cluster centre - and are 
isotropic or slightly tangential. Models for the gas velocity pro- 
files are required because the gas on average retains residual mo- 
tions even inside the virial radius, and this provides non-negligible 
anisotropic pressure support to the equilibrium of the system. For 
this reason we generalized the hydrostatic equilibrium equation by 
adding velocity terms analogous to those used in the Jeans equation 
for dark matter. This 'hydrodynamical equilibrium' equation is in 
fact a better description of the system and provides a more accurate 
estimator of the actual mass enclosed by any radius x. However, 
if gas velocity are neglected, a simple isothermal model fares bet- 
ter than a non-isothermal one: if our gas density profile is adopted, 
then the actual cluster mass is recovered to better than 15 per cent 
at all radii of interest. 

The dark matter and total density profiles in our SPH simula- 
tions appear to be more centrally concentrated than those coming 
from identical - but collisionless - simulations. In terms of NFW 
halo concentrations, these increase by an amount comparable to 
the baryonic fraction. Nevertheless, the addition of a hot gas com- 
ponent does not degrade (but even improves) the goodness of fit of 
a NFW profile. Therefore, the fact that our clusters are better fitted 
by the profile proposed in equation {3J is not due to the presence of 
the ICM. On the contrary we showed that it is due to the presence, 
in our sample, of clusters in various environmental situations, as 
opposed to having a sample of isolated systems. 

The shape of the gas density profile at small radii is at least 
partially explained by the gas expansion caused by the energy gain 
from dark matter during the collapse. Proper inclusion of the gas 
bulk energy shows that the final energy budget is equally shared be- 
tween the two species at radii r > 0.1W V , while at smaller radii the 
gas is hotter than the dark matter. This energy unbalance is proba- 
bly also the reason of the further global compression of the system 
compared to a pure collisionless collapse. It would be interesting to 
investigate whether this effect is large enough to produce observa- 
tional signatures. 

The analytical fits proposed in this work give an overall robust 
and accurate description of the average dark matter and hot gas 
properties of galaxy clusters. They have immediate applications in 
cosmology, as one can derive from them the radial dependence of 
many observables in the X-ray, SZ and lensing domains. 

Work is also in progress to extend this analysis to a sample 
of cluster simulations where the gas is heated by non-gravitational 
processes at high redshift, providing better agreement between the 
X-ray properties of simulated and observed clusters. 
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